
## Load packages 
library(readxl)
library(ggplot2) 
library(tidyverse)
library(ggtext)

## Read & sort data
R_import <- read_excel("[paste path]/R import comorb Asia.xls")

unique_countries <- unique(R_import$country)
unique_countries <- c(setdiff(unique_countries, "Global"), "Global")
R_import_ord <- R_import
R_import_ord$country <- factor(R_import$country, levels = rev(unique_countries))

R_import_ord <- R_import %>%
  arrange(region, Estimate) %>%  
  mutate(country = factor(country, levels = unique(country)))  

R_import_ord$type <- ifelse(R_import_ord$country == "GLOBAL", "GLOBAL", "Other")
R_import_ord$region <- as.factor(R_import_ord$region)

## Define labels for regions
region_labels <- c(
  "1" = "Europe & North America",
  "2" = "Central Asia",
  "3" = "Middle East & North Africa",
  "4" = "Latin America & Caribbean",
  "5" = "Pacific island nations",
  "6" = "South, East & Southeast Asia",
  "7" = "Sub-Saharan Africa",
  "8" = "Global"
)

## Create the plot        
forest <- ggplot(R_import_ord, aes(x = Estimate, y = fct_rev(country))) +
  geom_errorbarh(aes(xmin = LCL, xmax = UCL, color = region), height = 0.2) +
  geom_point(aes(color = region, size = type, shape = type)) +
  labs(x = "Share of eligibility due to comorbidity (%)", color = "World region") +
  theme_classic() +
  scale_color_hue(labels = region_labels) +
  scale_shape_manual(values = c("GLOBAL" = 18, "Other" = 16), guide = "none") +       # Custom shapes
  scale_size_manual(values = c("GLOBAL" = 4, "Other" = 2.5), guide = "none") +          # Custom sizes
  theme(
    legend.position = c(0.95, 0.95),  
    legend.justification = c(1, 1), 
    legend.box.just = "right",      
    legend.text = element_text(size = 10),      
    legend.title = element_text(size = 11, face="bold"), 
    legend.background = element_rect(fill = "white", color = "black"),  
    panel.grid.major.x = element_line(color = "lightgrey", linewidth = 0.2),
    axis.title.y = element_blank(),
    axis.text.y = element_text(color = "black", size = 10),
    axis.text.x = element_text(color = "black", size = 11),
    axis.title.x = element_text(color = "black", size = 11, face = "bold")) +
  scale_x_continuous(limits = c(0, 100), breaks = seq(0, 100, by = 10), expand = expansion(add = c(0, 2)))

ggsave("Comorbidity.pdf", plot = forest, width = 9, height = 13)

